
### This code generates the bootstrap simulations bs_simulation.X.RData
### Code is run in parallel on Tukey
### The output is later aggregated and plotted in binIRT simulation.R

rm(list = ls())
library(emIRT)
chain.id <- as.integer(Sys.getenv("SLURM_ARRAY_TASK_ID"))
set.seed(123L * chain.id)

## Set number of trials, and bootstrap iterations per trial, along with containers
N.trials <- 25	#1000 total, divided by processes
N.bs <- 100		#100 always

## Read original data set, run estimates
#h112 <- readKH("ftp://voteview.com/hou112kh.ord")
load("h112.rda")
rc <- convertRC(h112)
p <- makePriors(rc$n, rc$m, 1)
s <- getStarts(rc$n, rc$m, 1)
p1 <- p
p$x$sigma <- matrix(10,1,1)
rc.copy <- rc
NA.vec <- which(rc$votes %in% c(0,9))
h112.out <- binIRT(.rc = rc, .starts = s, .priors = p,
      .control = {list(threads = 10, verbose = FALSE, thresh = 1e-6) })

## Extract vote probabilities
voteprob.true <- pnorm(cbind(1,h112.out$means$x) %*% t(h112.out$means$beta))
idealpt.true <- (h112.out$means$x - mean(h112.out$means$x))/sd(h112.out$means$x)
N <- nrow(voteprob.true)
J <- ncol(voteprob.true)

## Initialize containers
bse.obs <- matrix(NA, N, N.trials)
idealpts.obs<- matrix(NA, N, N.trials)
coverage <- matrix(0, N, N.trials)
idealpts.bs<- matrix(NA, N, N.bs)
rownames(bse.obs) <- rownames(rc$votes)
rownames(idealpts.obs) <- rownames(rc$votes)
rownames(coverage) <- rownames(rc$votes)
rownames(idealpts.bs) <- rownames(rc$votes)

for(i in 1:N.trials){

## Estimate one "observed" realization of the data
votemat <- 1*(voteprob.true > matrix(runif(N*J), N, J))
votemat[votemat[[i]]==0] <- -1
votemat[NA.vec] <- 0
rc.copy$votes <- votemat
obs.out <- binIRT(.rc = rc.copy, .starts = s, .priors = p,
      .control = {list(threads = 1, verbose = FALSE, thresh = 1e-6) })

## Calculate vote probabilities for observed realization, store ideal points
voteprob.obs <- pnorm(cbind(1,obs.out$means$x) %*% t(obs.out$means$beta))
idealpts.obs[,i] <- (obs.out$means$x - mean(obs.out$means$x))/sd(obs.out$means$x)

	## Bootstrap the observed realization of data
	for(j in 1:N.bs){

	votemat.bs <- 1*(voteprob.obs > matrix(runif(N*J), N, J))
	votemat.bs[votemat.bs==0] <- -1
	votemat.bs[NA.vec] <- 0
	rc.copy$votes <- votemat.bs
	bs.out <- binIRT(.rc = rc.copy, .starts = s, .priors = p,
	      .control = {list(threads = 10, verbose = FALSE, thresh = 1e-6) })
	idealpts.bs[,j] <- (bs.out$means$x - mean(bs.out$means$x))/sd(bs.out$means$x)
	} #for(j in 1:N.bs)

bse.obs[,i] <- apply(idealpts.bs, 1, sd)

cat("\n\t\t ====================== Trial", i, "complete ==========================")
flush.console()

} #for(i in 1:Ntrials))


upper.95 <- idealpts.obs + 1.96*bse.obs
lower.95 <- idealpts.obs - 1.96*bse.obs
for(i in 1:nrow(coverage)){
  for(j in 1:ncol(coverage)){
	if( (upper.95[i,j] > idealpt.true[i]) & (lower.95[i,j] < idealpt.true[i])) coverage[i,j] <- 1
  }
}

votes <- matrix(1,N,J)
votes[NA.vec] <- 0
pc.vote <- apply(votes,1,sum)/J

save(idealpts.obs, bse.obs, idealpt.true, pc.vote, file = paste("bs_simulation", chain.id, "RData", sep = "."))

